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This paper studies the small sample behavior and convergence properties of 
a number of confidence interval estimators (CIE's) for the mean p of a 
stationary process, X^,...,X n . These CIE's are typically of the form 

Pr ^ 6 \ s 4 1 - «’ <*> 

where X = E.X./n, t . Y is the Y-quantile of the t-distribution with d degrees 
n ii ci p * 

of freedom (d.o.f.), and V estimates = nVar(X fi ) (or the variance parameter 

2 2 2 2 
a = lim a ). A "good" estimator for o (or o ) is the cornerstone of a 
n-«D n' b n' 

valid CIE for p. Many estimators have been studied in the context of discrete- 
event simulation: nonoverlapping batch means (NOBM) (Conway 1963, Schmeiser 
1982, Kang 1984); independent replications; overlapping batch means (OBM) 

(Meketon 1980, Meketon and Schmeiser 1984); standardized time series (STS) 

(Schruben 1983, Goldsman 1984, Glynn and Iglehart 1990); spectral theory 
(Fishman 1973,1978, Heidelberger and Welch 1981,1983); ARMA time series 
modeling (Fishman 1973,1978, Schriber and Andrews 1984); and regeneration 
(Crane and Iglehart 1975, Crane and Lemoine 1977, Fishman 1978). 

There is considerable work which compares the various CIE methodologies. 

The Monte Carlo (MC) work mainly deals with small sample CIE performance; see, 
e.g.. Law and Kelton (1984) and Goldsman, Kang, and Sargent (1986). Analytical 
results are almost all asymptotic: Goldsman and Schruben (1984), Goldsman and 
Meketon (1986), Damerdji (1987), Glynn and Iglehart (1990), and Schmeiser and 
Song (1989) all compare various combinations of the CIE's. 

In this paper we study finite sample behavior of CIE's from NOBM, OBM, and 
STS. Section 1 gives background material. Section 2 reports on statistical 
properties of various variance estimators. Section 3 presents analytical 
results for some special cases and then summarizes a MC study of the CIE's. In 
Part I of our MC work, we break the n observations into b batches, and we 



- 2 - 



observe what happens as the batch size grows. For "small" b, NOBM performs the 
best with respect to CIE coverage. For "large" b, both NOBM and OBM fare the 
best. However, an STS method produces CIE's with smaller expected lengths. 
Another comparison is carried out in Part II of our MC work, where we fix the 
d.o.f. d and observe what happens as n grows; here, the NOBM, OBM, and STS 
"combined" CIE's perform similarly. Section 4 discusses our findings, and 
Section 5 summarizes. We conclude that the bias of V is the most important 
factor in determining a CIE's validity; a secondary role is played by the 
marginal distribution of the X^'s. We also find that a CIE having superior 
large sample properties may have relatively poor small sample performance. We 
offer practical and research recommendations. 

1. BACKGROUND 

We review the CIE's and stochastic processes under study. We assume the 

stochastic processes satisfy certain moment and mixing conditions, as described 

2 2 

in the cited papers. We henceforth use the notation Nor(p,r ), X (d), X(d), 
Exp(X), and t(d) to represent the normal, chi-square, chi, exponential, and t 
distributions, each with the appropriate parameters. 

1 . 1 Nonoverlapping Batch Means 

Suppose we divide Xj,...,X r into b > 1 adjacent, nonoverlapping batches of 

size m (assume n = bm). The i-th batch mean, i = l,...,b, is X. = 
v ' i,m 

1° implementing the method of NOBM, we assume the X^ m 's are 

2 2 
approximately i.i.d. Nor(p,a /m) . The NOBM estimator for a is 

V N = m ^ =1 ( X. >m - X n ] 2 /(b-l) l a 2 X 2 ( b-l)/(b-l), 

2 ) 

where " ■+ " denotes convergence in distribution as m ■>«. The NOBM CIE for p 
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is given by (l) with d = b— 1 and V = V^. 

1.2 Overlapping Batch Means 

Define the i-th overlapping batch mean , i = 1, — ,n-m+l, by X(i,m) = 
y m i X., ./m. The OBM estimator for a 2 is 

"J=0 i+j' 

V Q = nm ^~“ +1 [ X(i,m) - X n ] 2 /[(n-m+l)(n-m)] , 

and is almost identical to Bartlett's spectral estimator (see Priestley 1982). 

The OBM CIE for p is given by (1) with V = V^; its validity depends on being 
2 2 

approximately a X (d)/d. Meketon and Schmeiser (1984) take d = 1.5* (b-1), where 
b = n/m. Based on MC experimentation, Schmeiser (1986) recommends d = 
1.5*(b-l)[l+(b-l) (*5+*6b)j. we ghall use this value in our MC work. 



1.3 Standardized Time Series 

Suppose X r ...,X n is divided into b ;> 1 adjacent, nonoverlapping batches 
of size m. For i = l,...,b, let A. = [(m+l)/2 - Schruben 



(1983) assumes the A. 's are approximately i.i.d. normal, and proposes the area 

2 

and combined estimators for a : 



V A^ 



12 



(m 2 -m)b 



, o ~ 2 v 2r, \ (b-l)V M + bV A 2 v 2^ ou .n 

yb j 2 5 £_xj a 5 —4J-, — i . 

£*i=l i b C 



2b-l 



2b-l 



CIE's for p are formed by substituting the appropriate V and d in (1). Schruben 



also derives the so-called "maximum" estimator for a (see Subsection 4.2). 



1.4 Some Time Series Processes of Interest 

fe will have occasion to use the following ARMA-type processes. 



MA(l) : X. = c. + 0£. „, 
l l l-l 

MA'(l): X. = £. + Be. . 

l l l-l 

AR(l) : X. = tpX. . + £., 

l T l-l i 



where c. ~ i.i.d. Nor(0,l) and -1 < 0 < 1. 
where ~ i.i.d. Exp(l) and -1 < 0 < 1. 



~ i.i.d. Nor(0,l-<p ), — 1 < (p < 1 . 



where £. 
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EAR(l) : X. = < 



<px 



i-1 



w.p. ip 



• vX i-i + e i w * p * 1_<p 



where ~ i.i.d. Exp(l) and 
0 £ ip < 1 (see Lewis 1980). 



We also consider other stationary normal processes as well as the waiting time 



(delay) processes of M/M/1 and E^/M/l queues. 



1.5 CIE Performance Criteria 

Denote the NOBM, OBM, area, and combined CIE's by CIE^, CIEq, CIE a> and 
CIEp, resp. The half-length of a generic CIE is H = t^ ^ _ a / 2 (V/ n )^- We use 
the following CIE performance measures: coverage (Pr$p E i H|), E[H], and 
Var(H). Among CIE's which achieve coverage 1 - a, we prefer that with the 
smallest E[H], and then that with the smallest Var(H). 



2. PROPERTIES OF VARIOUS VARIANCE ESTIMATORS 

We give some results on the bias and variance of the NOBM, OBM, and STS 
estimators, and on the asymptotic performance of the corresponding CIE’s. 

2.1 Bias of the Estimators 

The bias of V as an estimator for a is Bias(V) = o - E[V]. Goldsman and 

Meketon (1986) show that lim mlim, Bias(V„) = lim mlim, Bias(V_) = 
v 7 m->© b-»© N m->® b->© 0 

lim mlim. Bias(V„)/2 = lim mBias(V.)/3. All of these estimators are 
dh® b-»© C m-»® A 

asymptotically unbiased as m ■+ ®, but for small m, these estimators can be 
quite biased. 

p 

Example 1 : For the AR(l) and EAR(l), the Appendix gives o = (l+<p)/(l-<p) and 

E [V ] = a 2 - 2( P( b+1 > + 2<Pb(<p m -b~V b ) £ Z _ 2y 

^ mb(l-ip) 2 m(b-l) (l-<p) 2 m(l-(p) 2 



for large m and b, (2) 
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f - 



2bip 



m(b-l) (1-tp)' 



mb 



i - *" + 



] t «,dV)(i 

-I m(b-l)(mb-: 



m N/4 mb-m+l s 

± L 



2 

= a - 



2(p 



m(l-(p)' 



for large m and b. 



(b-l)(mb-m+l) (l-tp)' 



E [V 1 = , g 2 6j_ (4) 

' — ( 1 —in") L J m f 1 — a 



(3) 



E[V_] = <T + 



. 2 
= a - 



(m -m)(l-<p) 

2(£ 

(2b-l)m(l-<p) 
4<p 



(1-<P) 

/ _ 4b + 1 ~ < P mb - 3b(m+2)<p m 12b<p[l-(m+l-nnp)<p m ] 1 

1 b “- 1 (m 2 -l)(l-,) 2 J 



m(l-<p) 
for large m. 



(5) 



for large m and b. 



m(l-<py 

If m and b are large, the bias results anticipated by Goldsman and Meketon are 
attained. Table 1 contains exact E[V]'s for b = 2 and 16 and various m. For 
small b, Bias(V^) < Bias(V^) < Bias(V^) < Bias(V A ); for large b, Bias(V^) = 
BiasCV^) < Bias(V^) < Bias(V^). The biases become negligible as m grows. ^ 



2 2 

Example 2 : For the MA(1) and MA'(l), the Appendix gives a = (l+9) and 



E[V n ] = ct 2 - 20(b+l)/mb = a 2 - 29/m for large b. (6) 

e[v„] - <? - _ sfsfi] 4 _ 2e/m for laree b - (7) 

E[V a ] = a 2 - 60/m. (8) 

E[Vg] = a 2 - (4b+2)0/mb = a 2 - 49/m for large b. (9) 



The conclusions from Example 1 again hold. ^ 

Although Bias(V) is interesting in its own right, the bias directly 
affects CIE performance. Consider the unrealistic case in which m is fixed and 
b oo. For the estimators studied here, one can show that as b •» co, V E[V] 
w.p.l, and so T = (X n -p)/(V/n) X ® Nor(0,cr 2 /E[V]) . Falsely assuming T ~ 

Nor(0,l) as b o yields incorrect CIE coverage. 
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Pr||i 6 X n ± H* h. 2*[z 1 _ a/2 (E[\]/o Z ) y ’'] ~ 1 > ( 10 ) 

^ 2 

where $(•) and z y are the Nor(O.l) c.d.f. and Y-quantile. If E[V]/a < [>] 1, 
then coverage < [>] 1 - a. (Coverage is quite sensitive to decreases in 
E[V]/a .) So the less bias the better. As b ■+ ®, CIE^ and CIEq ^ en( ^ 
achieve the desired coverage more quickly with respect to m than do CIE A and 
CIE^,; see Sargent, Kang, and Goldsman (1987) (S-K-G). 

2.2 Variance of the Estimators 

Goldsman and Meketon report that as m and b become large, b-Var(V^) •* 2 a , 
b*Var(V Q ) ■+ 4a 4 /3, b-Var(V A ) ■+ 2a 4 , and b*Var(V c ) -► a 4 (cf. Damerdji 1987). 

For i.i.d. X^, — .X^, Kang and Goldsman (1990) find Var(V^) and Var(V A ) 
exactly, as do Song and Schmeiser (1989) for Var(V^). Exact results for other 
processes and for Var(V^) are tedious to derive. Of course, one can also 
calculate the mean squared errors of the V's (cf. Schmeiser and Song 1989). 



2.3 Asymptotic Properties of the CIE's 

Schmeiser (1982) and Goldsman and Schruben (1984) note that as m •+ ®, 



(mb)*H 5 <* d>1 _ a/2 *(d)/Md (11) 

for CIE^, CIE^, and CIE^. An analogous approximate result holds for CIEq. 

Under (ll), the CIE's achieve coverage 1 - a as m -► ®. Further, if r(«) is the 
gamma function, then 



("b)*S[H] - < rt d,i-«/2 (2/d),! r( r(d/2{ 2) • and 



mbVar(H) * { 1 



2 f r((d+i) 
d [ r(d/2 




( 12 ) 



(13) 



The right sides of (12) and (13) decrease in d and, hence, in b. So for large 
m with fixed b, E[H N ] > E[H A ] > E[H Q ] > E[H C ] and Var(H N ) > Var(H A ) > Var(H Q ) > 
Var(H^), the subscripts having the obvious meanings. Goldsman and Schruben 
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(1984) and Meketon and Schmeiser (1984) let b -► ® in (12) and (13) to get 



lim 



E[H] 






= 1, H = H 0 , H a> H c , 



and 



H ”W 

m-+® 



Var(H) 

Var(H N ) 



' 1, h = h a 

- 2/3, H = H q . 
.1/2, H = H C 



( 14 ) 



Thus, as b also becomes large, all of the CIE's have about the same E[H]'s, but 
CIEp has the smallest Var(H). 



3. FINITE SAMPLE CONFIDENCE INTERVAL ESTIMATION 



Small sample analysis of CIE's is difficult. We present a few exact 
results, but most of the section is devoted to a MC study. 



3. 1 Some Analytical Results 



Example 3 : Suppose X^,...,X ~ i.i.d. Nor(p,T^). Then ~ T^X^(b-l)/(b-l) , 

V A ~ T^X^(b)/b, and ~ r^X^(2b-l)/(2b-l) . Further, is independent of V^, 
V A> and (see Appendix). So (1) holds exactly for CIE^, CIE A> and CIE^. We 
could not obtain such results for CIEq or for nonnormal i.i.d. processes. ^ 



Example 4 : We can derive exact results for CIE A when b = 1 (n = m) and X^,...,X 

is stationary normal . Then A. ~ Nor(0,E[A^l) , V. ~ E[V. 1x^(1), and X is 

l l A A n 

normal and independent of V A (see Appendix). So (X^-p) (nc/V A ) * ~ t(l), where 

^ 2 

c = E[V.]/ct . Hence, the coverage of CIE, is 2Pr$t(l)<;t 4 . /0 *Jc^-l = 

An A 1 , 1— OC/ (Z 

(2/n)Tan" 1 (t 1 ^_ a yg-Jc). (As in (10), the coverage is sensitive to decreases 
in c.) Similarly, E[H A ] = t^ ]_ a /2^E[V^]/nn) *. To iH us trate, suppose the 
X.'s are AR(1). Figure 1(a) uses Example 1 and (A-3) to plot coverage vs. 
log^n for 1 - Of = 0.90 and various <p. We see that for <p > [<] 0, the coverage 
is < [>] 1 - a. As |<p| approaches 0 or as n grows, Bias(V A ) decreases, and the 
coverage approaches 1 - a. Figure 1(b) has analogous plots of E[H A ] vs. 
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A 2 —y 

loggn. If <p = 0.0, then E[V^] = a =1, and so E[H^] decreases at rate n *. 

If <p = -0.9, then E[V^] decreases to cP = 1/19, and E[H^] decreases at rate n * 

(after initially decreasing faster). The <p = 0.9 plot for E [H^] increases and 

then decreases as n grows. This occurs since E[V^] increases to a =19 as 

n ■+ ®, while the 'Jn in E[H^]'s denominator becomes large. ^ 

Example 5 : We give exact results for CIE^ when b = 2 (n = 2m) and X^, ,X R is 

stationary normal. Then V M = mE.(X. -X )^/(b-l) = m(X. - X„ )^/2; so V M ~ 

E[V^]X ,i (l). Since X^ is normal and independent of for b = 2 (see Appendix), 
we have (X R -p) (nc '/V^) * ~ t(l), where c’ = E[V^]/a n . The coverage is 
(2A)Tan“ 1 (t 1>1 _ a/2 >lc') > and E[H N ] = t 1>1 _ a / 2 (E[V N ]/nm) ;f . If the X.'s are 
AR(1), then (A-3) yields o^, and (2) with b = 2 gives E[V^]; c', coverage, and 
E[H^] then follow. These performance measures behave as in Example 4. j j 

It is difficult to generalize the above CIE results to b > 2, since we 

2 

would then have correlated X random variables. Thus, we only gave exact 
results for simple cases. We resort to MC experimentation in the sequel. 

3.2 Design of the Monte Carlo Study 

Our goal was to assess CIE performance over a variety of stochastic 

processes and choices of number of observations n, batch size m, d.o.f. d, and 

desired coverage 1 - a. We simulated the following processes: AR(1) with 
<p = 0.0, +0. 1, ±0.5, +0.9; EAR(l) with <p = 0.0, 0.1, 0.5, 0.9; MA(l) and MA'(l) with 
0 = ±0.1, ±0.9; M/M/1 with traffic intensity p = 0.6, 0.8 (service rate = 1.0); 
and Eg/M/1 with p = 0.6 (service rate = 1.0). Each run was initialized from the 
appropriate steady state distribution. All uniform [normal] variates were 
generated from algorithm UNIF [TRPNRM] in Bratley, Fox, and Schrage (1987); 
exponentials used inversion. 
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For ease of exposition, we divided the study into two parts. In Part I we 
fixed 1 - a and b, and then charted CIE coverage as a function of m. Roughly 
speaking, we wanted to know which CIE first achieved acceptable coverage as m 
increased with fixed b. Further, which CIE's had the best E[H] and (to a 
lesser extent) Var(H)? For the stochastic processes discussed above, we 
conducted at least 1000 independent runs of 16384 observations; these runs were 
used to calculate CIE^, CIEq, CIE a> and CIE C and the resulting performance 
characteristics for all choices of m = 2 , k = 0,1, — ,10, b = 1,2,4,8,16, and 
1 - a = 0.80,0.90,0.95,0.99. 

In Part II, we set 1 — oe = 0.90 and d = 3 and 15, and then charted 

coverage as a function of n. For d = 3 [d = 15], CIE^ uses b = 4 [b = 16], 

CIE,. uses b = 2 [b = 11], CIE. uses b = 3 [b = 15], and CIE P uses b = 2 

[b = 8] . We conducted our experiments on a number of stochastic processes, 

each of which used 2000 runs of 16384 observations to calculate the CIE's for 

all n = 2^, k = 4, — ,14, and d = 3 and 15. Since d is fixed, ( 11)— ( 13) 

suggest that all of the CIE's will perform about the same for large m; however, 

we suspected that they would exhibit different small sample performance because 

the variance estimators incorporated in the CIE's operate under different 

assumptions. We first discuss the underlying assumptions for V M , V. , and V_ 

since these estimators require independence between batches (Vq does not). The 

estimator V*. [V.] assumes that the X. 's [A.'s] are i.i.d. normal; the 
N L A J i,m L i J 

combined estimator V^, must satisfy both assumptions. For fixed d and n, and 
use roughly half the batch size of V^; hence, the assumption of i.i.d. 

X. 's [A.'s] is harder to achieve for V N [V ] than for V p . On the other hand, 

1 y ID 1 IN A L 

V^'s assumption of normality of the X.. m 's is easier to satisfy than V A 's 



i,m 



assumption of normality of the A^'s which relies on a more restrictive 
functional central limit theorem. The normality question for vs. is not 
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as clear since uses half of Vq' s batch size. CIEq appeals to spectral 
theory to directly assume that Vq is o X (d)/d; this supposition is not true 
for nonnormal processes or for finite batch sizes. 

3.3 Results from Part I of the Monte Carlo Study 
3.3.1 Representative results 

We discuss typical results from Part I. Figures 2, 3, and 4 illustrate 
CIE performance when X^, — ,X fi are AR(1) (<p = 0.9), EAR(l) (<p = 0.9), and M/M/1 
(p = 0.8), resp. , and 1 - a = 0.90. The AR(l) and EAR(l) have the same 
covariance function, but the AR(1) has normal marginals while the EAR(l)'s are 
exponential. The M/M/l's joint distribution is more complicated. We only 
consider the cases b = 2 and 16 since CIE^, CIEq, and CIEq require b ;> 2, and 
since one can argue that b = 16 is "large" (at least for the "usual" choices 
of a). Each of Figures 2, 3, and 4 has four graphs: (a),(b) are for b = 2, and 

(c),(d) are for b = 16. In (a) and (c), we plot the sample coverage (CVG) vs. 

log^m; (b) and (d) plot the sample E[H] (EHL) vs. loggin. The standard error of 
any CVG is about [CVG* (1-CVG)/1000]^. 

Figure 2 is for the AR(l) with (p = 0.9. All CVG's are poor for small m, 
but approach 1 - a = 0.90 as m increases. For b = 2, the CVG of CIE^ is 
closest to 0.90; for b = 16, both CIE^ and CIEq yield the best CVG's. This 
makes sense since V^ (for b = 2 and 16) and Vq (for b = 16) are less biased 

than V^ and Vq (Example 1 and Table 1). A related consequence is that CIE^ and 

CIEq produce larger EHL's than those of CIEq. Note that the EHL's in Figure 2 
increase and then decrease as m increases - the same bias-related pattern as for 
the exact E[H^] in Figure 1. The EHL’s for b = 16 and m ;> 128 from Figure 2(d) 
are more or less in agreement with the limiting E[H]/E[H^] = 1 ratios in (14); 
this is one reason why we regard b = 16 as "large." 
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Table 2 has CVG's for the AR(1) with b = 16, 1 - a = 0.90, and various 

m and <p. These roughly agree with the b ■+ ® results from (10). We see that 

CVG < [>] 0.90 when <p > [<] 0.0. For any <p and fixed m. Table 2 shows that the 

CVG's of CIE^ and CIEq are closer to 1 - a than those of CIE A and CIE^. 

Results for the EAR(l) process with <p = 0.9 are found in Figure 3, whose 

plots bear a striking resemblance to those of the AR(l) in Figure 2. The only 

notable difference between the two figures is that, for fixed m, the EAR(l) has 

smaller CVG's (see Table 1). Since the AR(l) and EAR(l) have the same 

covariance structure, the EAR(l)’s poorer CVG's are probably due to its 

exponential marginals. It seems that the effect of Bias(V) on coverage is more 

significant than that of the marginals of the X^'s. 

Figure 4 concerns the M/M/1 process with p = 0.8. Again, the patterns in 

the figure are not much different than those of the AR(l) and EAR(l). The 

M/M/1 simply requires more observations to attain valid coverage. The positive 

2 

serial correlation of the M/M/1 causes the estimators for a to be biased too 
low; so poor coverage results for small m. 

3.3.2 Additional Part I results 

In S-K-G, we also give detailed discussions on (among other things): 

- CIE performance for the MA(l) as 0 varies, fe find that the CIE's 
qualitatively perform about the same as those for the AR(l). 

- The sample -JVar(H) (SHL) performance measure. For small m, the SHL's exhibit 
the same general behavior as the EHL's; as m and b become large, the SHL's 
behave as in (14). Even though the ratios from (14) are manifested, the 
differences between SHL's from competing methods are typically very small. 
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- The consequences of changing 1 - a. The CIE's qualitatively perform the same 
as we vary 1 - a. Coverage is sensitive to decreases in 1 - a. For example, 
consider CIE^ for the AR(l) with <p = 0.9, b = 2, and m = 16; for 1 - ot = 0.99, 
0.90, and 0.80, Example 5 gives coverages of 0.985, 0.852, and 0.713, resp. 

3.3.3 Summary of Part I 

For small m, improper CVG's were usually the rule. For small m and b, 

CIE^ has better CVG's than the other CIE's, while for small m and large b, both 
CIE^ and CIEq fared the best. For fixed m, high CVG was most often accompanied 
by high EHL and SHL. The CIE's performed as expected by asymptotic theory when 
m and/or b became large. For large m and small b, all achieved the desired 
coverage, and the EHL's and SHL's tended to decrease with increasing d.o.f., as 
per Subsection 2.3. For large m and b, the ratios from (14) took effect. 

3.4 Results from Part II of the Monte Carlo Study 

We considered the MA(l) (0 = -0.9), AR(l) (<p = 0.9), EAR(l) (<p = 0.9), and 

M/M/1 (p = 0.8), with 1 - a = 0.90 and d.o.f. d = 3 and 15. Table 3 gives CVG's 

as a function of the sample size n. For d = 3, CIE^, CIEq, and CIEq perform 

about the same in terms of CVG; CIE^ fares poorly for small n. For d = 15 and 

small n, CIE~ does a bit better than the others in terms of CVG; CIE. is not 
U A 

competitive. However, the performance of CIEq is "more variable" than the 
other CIE's over the range of stochastic processes, d, and n under study. For 
instance, the CVG of CIEq sometimes significantly overshoots 1 - a, especially 
for small d (though this is understandable since 0BM was designed for 
large b) . As n grows (for d = 3 or 15), it appears that CIE^, CIEq, and CIEq 
achieve CVG = 0.90 at about the same n. 
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4. DISCUSSION 

We first discuss the causes of improper CIE coverage; why are some CIE's 
better than others? We then consider the question of which CIE is "best?" 

4.1 Causes of Improper Coverage 

A CIE of the form in (l) will attain perfect coverage if its associated 
pivot T = - p)/(V/n)^ ~ t(d); this requires (i) 3^ ~ Nor(fi,a^/n) , (ii) V ~ 

a^X (d)/d, and (iii) independence of and V. Requirement (i) is satisfied if 
the marginal distribution of the X^'s is normal. If the X^’s are not symmetric, 
then T may be skewed for small n. But in most cases, a central limit theorem 
asserts that (i) approximately holds as n grows. We believe that (ii) is the 
key requirement. At a minimum, V must be approximately unbiased as an 
estimator of a (or a^) . In fact, since variance directly affects the CIE's 
length, we claim that Bias(V) is often the main cause of improper coverage (at 
least for small m) ; see below. Concerning (iii), Glynn (1982) and Kang and 
Goldsman (1990) demonstrate that a symmetry in coverage is directly related to 

ZZ A, 

dependence between X^ and V. However, Kang and Goldsman give examples which 
show that actual coverage is not necessarily affected by such dependence. We do 

“ -A 

not view dependence between X r and V as a direct cause of improper coverage. 

We first analyze the effect of Bias(V) on requirement (ii) by examining V = 
a^V/E[V] instead of V; note that E[V] = . To illustrate, we shall use the 

NOBM estimator on the AR(1) and EAR(l) processes with <p = 0.9. For the AR(l) 
with b = 2, Example 5 says that ~ E[V^]X^(1), and so ~ a^X^(l); for this 
case, correction for bias results in precisely the desired distribution. 
Empirical p.d.f.'s of and (based on 100000 independent runs) are 

plotted in Figure 5 for the AR(l) and EAR(l) with b = 8 and various m. For the 
AR(l), the sample p.d.f.'s of v ^/ CT n approach the X (7) p.d.f. as m increases; 
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^ 2 2 

the corrected V^/a^ is nearly (but not quite) X (7) (Figures 5(a) and 5(b)). 

The empirical p.d.f.'s for the EAR(l) exhibit similar behavior, except that the 

2 

convergence to the X (7) is slower (Figures 5(c) and 5(d)). Thus, for these 
examples, correction for bias mitigates the violation of (ii). Comparing the 
AR(l) results with the EAR(l)’s, we conclude that the effect of Bias(V) is 
greater than that of the marginal distribution. 

We will next show for the stochastic processes investigated that 
correcting for Bias(V) results in good CIE's. The notation CCVG in Table 1 is 
the exact (from Example 5) or sample coverage obtained for the AR(l) and EAR(l) 
models with <p = 0.9 and b = 2 and 16. For the AR(1) with b = 2 and any m, the 
corrected NOBM pivot (S r - p)/(V^/n)^ ~ t(l). So the CCVG's for CIE^ are 
exactly 1 - oc = 0.90; thus, for this example, Bias(V^) is the sole cause of 
improper coverage. The corresponding MC CCVG improvements for the EAR(l) with 
b = 2 are significant but not as good as those for the AR(l), this indicating 
that a secondary marginal distributional effect is present. This conclusion is 
illustrated yet again by Figure 6. (The EAR(l)'s empirical p.d.f.'s are 
somewhat skewed for small m as explained earlier.) The results from Table 1 
and Figures 5 and 6 suggest that, for small m, there are also tertiary 
contributors to improper coverage, perhaps inter-batch correlation. 

4.2 Which is the Best CIE? 

The question of which CIE is the "best" depends on the criteria being 
used. A CIE is first judged by the validity of its coverage. But the MC work 
showed that a CIE with good coverage might produce relatively wide half- 
lengths; so the E[H] and Var(H) measures can not be ignored. Since coverage, 
E[H], and Var(H) are functions of the stochastic process, the CIE in use, the 
number of batches b, the batch size m, and the level a, we can see that the 
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determination of the "best" CIE is not straightforward. 

Results (12) and (13) say that E[H] and Var(H) decrease in the d.o.f. d as 
m becomes large with fixed b; so the more d.o.f. the better (although Schmeiser 
1982 finds that there is little to be gained by taking d > 30). Nevertheless, 
it would be incorrect to conclude that CIE^ (which has the largest d.o.f. of 
those CIE's under study) is always the best. The Part I MC results showed that 
for small b, CIE^ required the smallest value of m to achieve valid coverage 
(for large b, CIE^ and CIEq required the smallest m) ; if coverage were the only 
criterion for CIE comparison, CIE^ would be declared the best - not CIE^,. This 
shows that large sample superiority does not necessarily extend to the small 
sample case. Indeed, we did not include the STS "NOBM+maximum" CIE from 
Schruben (1983), which has 4b-l d.o.f. (and hence superior asymptotic 
properties), since it exhibited poor small sample performance compared to the 
other CIE’s (including CIE^). For instance, for the AR(l) with <p = 0.9, 1 - a 
0.90, b = 2, and m = 16, 64, and 256, the NOBM+max CIE attained CVG’s (based on 
1000 independent simulation runs) of 0.397, 0.606, and 0.776, resp.; CIE^ 
achieved CVG's of 0.684, 0.849, and 0.895, resp. (Table 1). 

Another basis for comparison among CIE’s is to determine which requires 
the smallest sample size n to achieve valid coverage for some fixed d.o.f. d. 
This was the aim of Part II of the MC study, where CIE^, CIEq, and CIE^ fared 
more or less the same; CIE. was not competitive. So there was no clear winner 
using the criterion of coverage under fixed d.o.f. 

We can still make some recommendations (for fixed sample size 
procedures). All of the CIE's studied here are easy to use. Batch means is 
the simplest method to understand. All are asymptotically valid as the batch 
size m grows; but it "never hurts too much" to use CIE^ (in comparison to the 
other CIE’s) in case m is not "large enough." The price to be paid when m is 
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small is that E[H^] and Var(H^) are larger than their competitors, particularly 
when b is also small. If the user is somehow confident that the batch sizes 
are large enough to achieve valid coverage, then the user should fix the d.o.f. 
(perhaps between 15 and 30 for the "usual" a values) and select from among 
CIE„, CIEq, and CIE^. However, one of the most difficult open questions in 
simulation output analysis is the determination of "sufficient batch size" (cf. 
Fishman 1978 and Schmeiser and Song 1989). 

5. SUMMARY AND CONCLUSIONS 

We studied the behavior of different CIE's with special emphasis placed on 
small sample size performance for various stationary stochastic processes. If 
the achieved coverages are at the desired levels, and if the total number of 
observations n is fixed, the ranking of the CIE's with respect to E[H] and 
Var(H) is determined by the d.o.f. each CIE has. In this case, the CIE with 
the largest d.o.f. has the smallest mean and most stable confidence interval 
length. 

Perhaps our most important finding was that, in small sample settings, a 
CIE with more d.o.f. may not actually be "better" than a competing CIE; some 
CIE's may require more observations than others before the asymptotics 
necessary for CIE validity manifest themselves. Quite often, the CIE's with 
the highest d.o.f. 's performed the most poorly in terms of coverage! 

The bias of V as an estimator of a (or a ) plays a significant role in 
CIE performance - the less bias the better. For instance, when m and b are 
fixed, the relative performance of the CIE's with respect to coverage is 
directly related to Bias(V). A secondary factor in CIE performance concerns 
the underlying marginal distribution of the X^’s. Further, with fixed m (and 
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even more so with fixed n) , coverage often deteriorates as b or a increase; 
this is partially due to the fact that, in these cases, t^ j_ a /2 decreases. 

Which CIE should one use in practice? If the sample size n is "large 
enough," we could probably argue successfully for CIE^, CIEq, or CIE^, with 
common d.o.f., 15 < d < 30. In comparison to the other methods, CIE^ is 
probably the "safest" small sample method. There are several interesting 
research lines. We would like to see more emphasis on small sample results, 
including sequential procedures (which were not investigated here). Another 
question concerns the fact that for fixed d.o.f., CIE^, CIEq, and CIE^, 
achieve approximately valid coverage for about the same sample size. 

Further, a good batch size estimation procedure would be of tremendous 
import . 
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APPENDIX 



Proof of (2)-(5) : We will use the following facts: 





(A-2) 



For the AR(l), and so 



2 = l+<e _ 2(p(l-(p n ) 

n ^ n(l-T) 2 




(A-3) 



and 
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Since 



V E ^N ] = ^j=l E[5 j,m ] " bE[ ^n ] = b ^[^ ] ' E[X n^ = b t Var (X m )-Var(X n )] , (A-4) 



result (2) follows by (A-3) and simplification (see Moran 1975). 



We also have E[ v 0 ] = ^ +1 E[(X(i,m) - J/] 



// 



r-m-1 r n 



:":f‘ E[x(i,m)x n ]. 


(A-5) 


e[x . , .x. ] = y? , v. . 

L i+j k J ^j=i ^k=l j-k 


(A-6) 



ym+i-1 ry j j-k + yn k-j-1 = ym+i-1 1 + y - <p J - y n J + 1 

“j=i ^**k=l v ^*k=j+l v * “j=i l-tp 



= m- 



(A-7) 



l+(p mw i . n-m+2 -i\ u . x2 

i-pj " (1-<P )(<P + <P <P )/(l-<p) • 

We obtain (3) by substituting (A-3) and (A-7) into (A-5) . ^ 

(4) follows from (A. 5-15) of Goldsman (1984); (5) follows from (2) and (4). ^ 



2 

Proof of (6)-(9) : The MA(l) has covariance function = 1 + 9 , Y + 1 = 6, and 
Y^ = 0, otherwise. By (A-2) , 

a 2 = (1+9) 2 - 29/n and a 2 = lim^^a 2 = (1+9) 2 . (A-8) 

Result (6) then follows from (A-4) and (A-8). ^ 



yn 

^k=l 



j-k 



’ (1+9) 2 

. (i+e) 2 -e 



if 1 < j < n 



if j = 1 or n 



Note that 
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So by (A-6) , 



mnE[X(i,m)X n ] = 



m(l+8)^ if 1 < i < n-m+1 

m(l+8)^-0 if j = 1 or n-m+1 



(A-9) 



We obtain (7) by substituting (A-8) and (A-9) into (A-5) . ^ 



(8) follows from (A. 5-23) of Goldsman (1984); (9) follows from (6) and (8). jj 

Proof of independence in Examples 3, 4. and 5 : Let X' = (X^,...,X n ), where the 

X^'s are stationary normal with covariance matrix E. Suppose G is an n x n 

symmetric matrix and V = X'GX. Problem 1.22 of Muirhead (1982) says that V and 

are independent iff l'ZG = O', where 1' [O'] is a lxn vector of l's [0's]. 

Song and Schmeiser (1989) note that V.. = X'GuX and V. = X'G.X, where 
& v/ N ~ N~ A A~ 




and G. = (h.h.) for b = 1, 
A i J 



where J is an n/2xn/2 matrix consisting of l's, and h. = (n+l)/2 - i. 



i 



i = 1 



> ♦ ♦ * y n * 



Then G^ and G^ both meet Muirhead's condition for any Z. ^ 
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Table 1 - Some Results for the AR(1) and EAR(l) Models with <p = 0.9. All 
entries for E[V] are exact. Results for CVG and CCVG are^exactC*) or are based 
on 2000 independent simulation runs. For these models, a = 19. 



E[V] AR(1) EAR(l) 

(or oj) b = 2 b = 16 b = 2 b = 16 





b=2 


b=16 


CVG 


CCVG 


CVG 


CCVG 


CVG 


CCVG 


CVG 


CCVG 


m = 4 






















NOBM 


0.86 


2.68 


0.744* 


0.900* 


0.479 


0.866 


0.625 


0.798 


0.465 


0.810 


OBM 


0.58 


2.64 


0.453 


0.840 


0.466 


0.859 


0.337 


0.673 


0.446 


0.804 


Area 


0.31 


0.31 


0.418 


0.907 


0.183 


0.888 


0.296 


0.701 


0.164 


0.799 


Comb 


0.49 


1.46 


0.432 


0.878 


0.365 


0.862 


0.321 


0.683 


0.337 


0.806 


( 2 
(CT n 


6.19 


16.19) 


















m = 16 






















NOBM 


6.10 


9.27 


0.852* 


0.900* 


0.756 


0.896 


0.791 


0.850 


0.742 


0.872 


OBM 


4.22 


9.23 


0.673 


0.864 


0.750 


0.890 


0.609 


0.773 


0.733 


0.866 


Area 


2.85 


2.85 


0.684 


0.896 


0.491 


0.891 


0.615 


0.800 


0.480 


0.866 


Comb 


3.94 


5.96 


0.689 


0.884 


0.651 


0.893 


0.625 


0.795 


0.646 


0.868 


, 2 
(a 
v n 


13.57 


18.30) 


















m = 64 






















NOBM 


14.79 


16.02 


0.891* 


0.900* 


0.874 


0.904 


0.881 


0.889 


0.869 


0.892 


OBM 


12.84 


16.00 


0.841 


0.883 


0.862 


0.894 


0.810 


0.841 


0.867 


0.888 


Area 


11.29 


11.29 


0.849 


0.897 


0.801 


0.894 


0.822 


0.866 


0.803 


0.884 


Comb 


12.45 


13.58 


0.852 


0.895 


0.840 


0.898 


0.819 


0.858 


0.841 


0.888 


( 2 
(a n 


17.59 


18.82) 


















m = 256 




















NOBM 


17.95 


18.25 


0.898* 


0.900* 


0.899 


0.905 


0.900 


0.902 


0.886 


0.890 


OBM 


17.30 


18.25 


0.895 


0.904 


0.902 


0.905 


0.900 


0.909 


0.890 


0.897 


Area 


16.90 


16.90 


0.895 


0.902 


0.886 


0.903 


0.888 


0.895 


0.878 


0.896 


Comb 


17.25 


17.56 


0.890 


0.901 


0.897 


0.908 


0.896 


0.902 


0.880 


0.892 


( 2 
(a 
v n 


18.65 


18.96) 


















m = 1024 




















NOBM 


18.74 


18.81 


0.900* 


0.900* 


0.897 


0.899 


0.883 


0.883 


0.909 


0.910 


OBM 


18.56 


18.81 


0.905 


0.907 


0.897 


0.898 


0.915 


0.918 


0.902 


0.903 


Area 


18.47 


18.47 


0.901 


0.905 


0.893 


0.897 


0.893 


0.895 


0.898 


0.901 


Comb 


18.56 


18.64 


0.892 


0.896 


0.900 


0.902 


0.895 


0.897 


0.899 


0.903 



(a|j 18.91 18.99) 
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Table 2 - CVG Results for the AR(l) Process, b = 16, 1 - a = 0.90. 
All entries are based on at least 1000 independent simulation runs. 



<p 


-0.9 


-0.5 


II 

a 


NOBM 


0.962 


0.945 


OBM 


0.966 


0.949 


Area 


1.000 


0.974 


Comb 


0.998 


0.971 


m = 16 


NOBM 


0.951 


0.915 


OBM 


0.958 


0.922 


Area 


0.988 


0.942 


Comb 


0.981 


0.926 


m = 64 


NOBM 


0.918 


0.904 


OBM 


0.929 


0.916 


Area 


0.943 


0.901 


Comb 


0.931 


0.900 


m = 256 


NOBM 


0.915 


0.903 


OBM 


0.915 


0.911 


Area 


0.931 


0.914 


Comb 


0.921 


0.918 


m = 1024 


NOBM 


0.906 


0.901 


OBM 


0.903 


0.899 


Area 


0.907 


0.904 


Comb 


0.909 


0.906 



0.0 


0.5 


0.9 


0.913 


0.823 


0.479 


0.908 


0.816 


0.466 


0.898 


0.677 


0.183 


0.900 


0.755 


0.365 



0.905 


0.887 


0.756 


0.906 


0.895 


0.750 


0.908 


0.860 


0.491 


0.910 


0.880 


0.651 



0.903 


0.901 


0.874 


0.910 


0.908 


0.862 


0.885 


0.885 


0.801 


0.896 


0.894 


0.840 



0.904 


0.900 


0.899 


0.909 


0.909 


0.902 


0.918 


0.917 


0.886 


0.916 


0.915 


0.897 



0.899 


0.899 


0.897 


0.897 


0.897 


0.897 


0.906 


0.904 


0.893 


0.902 


0.900 


0.900 
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Table 3 - CVG Results for Variance Estimators with Common d.o.f. 
All entries are based on 2000 independent simulation runs. 



n = 


32 


64 


128 


256 


512 


1024 


2048 


4096 


8192 


AR(l), 

NOBM 


to = 0.9 
0.688 


. d = 3 
0.787 


0.848 


0.880 


0.889 


0.899 


0.896 


0.902 


0.905 


OBM 


0.692 


0.792 


0.834 


0.874 


0.896 


0.910 


0.903 


0.923 


0.920 


Area 


0.505 


0.674 


0.795 


0.858 


0.880 


0.900 


0.895 


0.889 


0.912 


Comb 


0.701 


0.799 


0.847 


0.886 


0.883 


0.903 


0.895 


0.905 


0.908 


EARfl) 

NOBM 


. <P = 0- 
0.622 


9. d = 
0.736 


3 

0.826 


0.862 


0.887 


0.900 


0.896 


0.903 


0.896 


OBM 


0.631 


0.725 


0.820 


0.864 


0.887 


0.906 


0.917 


0.912 


0.916 


Area 


0.426 


0.630 


0.771 


0.844 


0.868 


0.886 


0.893 


0.893 


0.893 


Comb 


0.636 


0.740 


0.822 


0.863 


0.881 


0.898 


0.894 


0.902 


0.896 


M/M/1. 

NOBM 


O = 0.8 
0.437 


. d = 3 
0.541 


0.632 


0.707 


0.768 


0.804 


0.836 


0.862 


0.877 


OBM 


0.432 


0.540 


0.627 


0.704 


0.752 


0.799 


0.842 


0.867 


0.886 


Area 


0.290 


0.420 


0.533 


0.624 


0.708 


0.769 


0.817 


0.835 


0.855 


Comb 


0.442 


0.548 


0.634 


0.701 


0.758 


0.806 


0.838 


0.860 


0.871 


MA(1), 

NOBM 


6 = -0. 
0.985 


9. d = 
0.981 


3 

0.986 


0.975 


0.953 


0.932 


0.931 


0.920 


0.911 


OBM 


1.000 


1.000 


1.000 


1.000 


0.998 


0.986 


0.971 


0.959 


0.944 


Area 


0.994 


0.993 


0.989 


0.980 


0.967 


0.955 


0.947 


0.933 


0.918 


Comb 


0.994 


0.992 


0.989 


0.976 


0.963 


0.932 


0.933 


0.919 


0.908 



Aid), 

NOBM 


<0 = 0.9 
0.381 


, d = 
0.485 


15 

0.642 


0.749 


0.830 


0.877 


0.890 


0.900 


0.908 


OBM 


0.377 


0.515 


0.684 


0.793 


0.847 


0.880 


0.893 


0.897 


0.910 


Area 


0.126 


0.189 


0.318 


0.513 


0.695 


0.809 


0.871 


0.891 


0.903 


Comb 


0.385 


0.491 


0.649 


0.755 


0.834 


0.876 


0.894 


0.903 


0.907 


EARfl) 

NOBM 


. <P = Q- ! 
0.324 


9. d = 
0.470 


15 

0.607 


0.749 


0.811 


0.876 


0.903 


0.897 


0.897 


OBM 


0.319 


0.505 


0.660 


0.786 


0.833 


0.877 


0.899 


0.891 


0.902 


Area 


0.076 


0.173 


0.303 


0.496 


0.670 


0.802 


0.868 


0.882 


0.887 


Comb 


0.326 


0.479 


0.615 


0.754 


0.814 


0.879 


0.893 


0.894 


0.898 


M/M/1. 

NOBM 


o = 0.8 
0.227 


, d = 
0.318 


15 

0.386 


0.490 


0.618 


0.701 


0.773 


0.819 


0.857 


OBM 


0.226 


0.342 


0.425 


0.525 


0.653 


0.725 


0.779 


0.820 


0.858 


Area 


0.055 


0.091 


0.144 


0.228 


0.369 


0.520 


0.647 


0.757 


0.806 


Comb 


0.229 


0.319 


0.392 


0.498 


0.624 


0.706 


0.778 


0.823 


0.850 


MAfl) . 
NOBM 


9 = -0.9 
1.000 


d = 
1.000 


15 

1.000 


1.000 


0.999 


0.993 


0.983 


0.961 


0.948 


OBM 


1.000 


1.000 


1.000 


1.000 


1.000 


0.997 


0.983 


0.962 


0.952 


Area 


1.000 


1.000 


1.000 


1.000 


1.000 


1.000 


0.997 


0.990 


0.978 


Comb 


1.000 


1.000 


1.000 


1.000 


0.998 


0.992 


0.980 


0.963 


0.950 
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Figure 1: AR(l) Process Performance Me 
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Figure 3: EAR(l) Process Performance Measure 
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